#include "Molmodel.h"
#include <iostream>
#include <fstream>
using namespace SimTK;

int main() {
try {
   // Initialize simbody objects
	CompoundSystem system;
	SimbodyMatterSubsystem matter(system);
	DecorationSubsystem decoration(system);
	DuMMForceFieldSubsystem forces(system);
	forces.loadAmber99Parameters();
  
	// Initialize molmodel objects
	RNA rna("ACCAAUUAGGUUCAAAUAAGUGGU");  //ACGU
	rna.assignBiotypes();
	system.adoptCompound(rna);

	

	//rna.setRNABondMobility(BondMobility::Default, ResidueInfo::Index(0),ResidueInfo::Index(1));
	// Leave second residue at default combination of 
	// Torsion and Rigid mobilities
	

	//rna.setRNABondMobility(BondMobility::Default, ResidueInfo::Index(2),ResidueInfo::Index(2));
   
	// Assign atoms to bodies
	system.modelCompounds();

	// Maintain temperature
	system.addEventHandler(new VelocityRescalingThermostat(
           system,  293.15, 0.1));
   // Show me a movie
	Visualizer viz(system);
	//viz.setMode(Visualizer::Mode(2));
	system.addEventReporter( new Visualizer::Reporter(viz,.01));
  
	// Instantiate simbody model and get default state
	State state = system.realizeTopology();
   
	// Relax the structure before dynamics run
	LocalEnergyMinimizer::minimizeEnergy(system, state, 15.0);
 
	// Prepare for molecular dynamics
	VerletIntegrator integrator(system);
	integrator.setAccuracy(0.1);
	TimeStepper timeStepper(system, integrator);
	timeStepper.initialize(state);

	// Start simulation
	timeStepper.stepTo(500);
	return 0;
} 

catch(const std::exception& e) {
    std::cerr << "ERROR: " << e.what() << std::endl;
    return 1;
}
catch(...) {
	std::cerr << "ERROR: An unknown exception was raised"
              << std::endl;
    return 1;
}
}